Efficient simulation of one-dimensional quantum many-body systems 
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We present a numerical method to simulate the time evolution, according to a Hamiltonian made 
of local interactions, of quantum spin chains and systems alike. The efficiency of the scheme depends 
on the amount of the entanglement involved in the simulated evolution. Numerical analysis indicate 
that this method can be used, for instance, to efficiently compute time-dependent properties of low- 
energy dynamics of sufficiently regular but otherwise arbitrary one- dimensional quantum many-body 
systems. 
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Since the early days of Quantum Theory, progress in 
understanding the physics of quantum many-body sys- 
tems has been hindered by a serious, well-known compu- 
tational obstacle. The number of parameters required to 
describe an arbitrary state of n quantum systems grows 
exponentially with n, a fact that renders the simulation 
of generic quantum many-body dynamics intractable. 

A number of techniques have been developed to ana- 
lyze specific quantum many-body problems. Exactly and 
quasi-exactly solvable models [l| offer valuable insight 
in particular cases, and their solutions can be used as 
the basis for perturbative studies. With quantum Monte 
Carlo (QMC) calculations @, properties of the ground 
state of a large class of many-body Hamiltonians can be 
evaluated. In one dimensional settings, including quan- 
tum spin chains, ground-state expectation values for local 
observables, such as the energy and two-point correlation 
functions, can be computed with extraordinary accuracy 
with the density matrix renormalization group (DMRG) 
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The importance and ingenuity of these methods can- 
not be understated. However, solvable models and their 
extensions using perturbation theory apply to a very re- 
stricted class of physical systems, whereas QMC and 
DMRG calculations produce reliable outcomes mainly 
only for static properties of certain ground states. In par- 
ticular, the efficient simulation of time-dependent prop- 
erties remains an open problem in most non-perturbative 
cases. In addition to severely limiting our understanding 
of quantum collective phenomena, such as high- Tc super- 
conductivity and quantum phase transitions, the inabil- 
ity to efficiently simulate quantum dynamics has deep 
implications in Quantum Information Science Q. It is 
also a practical obtacle for the development of technol- 
ogy based on engineered quantum systems. 

In this paper we describe a numerical scheme for sim- 
ulating certain quantum many-body dynamics. We ex- 
plain how to efficiently simulate Hamiltonian evolutions 
in one-dimensional arrays of quantum systems, such as 
quantum spin chains, with a computational cost that 
grows linearly in the number n of systems. The key idea 
is to exploit the fact that in one spatial dimension, low 
energy quantum dynamics are often only slightly entan- 
gled. We employ a technique developed in the context 



of quantum computation [j| — thereby illustrating how 
tools from Quantum Information Science may find ap- 
plications in other areas of Quantum Physics (see also 
ESQ)- The present numerical method can be used, 
for instance, to determine the time-dependent vacuum 
expectation value 
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where O is a local Heisenberg operator. In some cases, 
and for specific choices of O, the Fourier transform 
S(k,oj) of is accessible experimentally, e.g. through 
neutron scattering for phonon dynamics in a solid. This 
allows for a direct comparison of experimental and sim- 
ulated data. 

We consider a one-dimensional array of quantum sys- 
tems, labeled by index I, I G {l,-- - ,n}, each one de- 
scribed by a local Hilbcrt space Hd of finite dimension 
d. We assume a collective evolution of the n systems ac- 
cording to a (possibly time-dependent) Hamiltonian H n 
made of local interactions, i.e. H n is a sum of terms each 
one involving at most k systems, for a small k indepen- 
dent of n. Physically reasonable Hamiltonians are of this 
form. Given a pure state \$>) 6 TLd^ n of the n systems, 
the entanglement between a block A containing the m 
first systems, I € {1, • • • , m}, and a block B containing 
the n—m remaining systems can be characterized by the 
rank Xa of the reduced density matrix pa for block A, 

X A = rank(p A ), Pa = tr B (|#)(#|). (2) 
Central in the present discussion is parameter X, 



X = maxX^, 



(3) 



the maximum, over the size m G {1, • • • ,n — 1} of block 
A, of the entanglement between blocks A and B [jj. 

The essential requirement for the simulation to be ef- 
ficient is that the entanglement X (in practice a related, 
effective parameter X € ) remains "small" during the sim- 
ulated dynamics, in a sense to be specified later on. Nu- 
merical tests indicate that this condition is usually met, 
for instance, when the Hamiltonian H n contains only 
short-range interactions and, in particular, in dynamics 
involving the ground state of H n and its low-energy exci- 
tations. Below we assume that H n contains precisely only 
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short-range interactions in an open chain, since these ad- 
ditional, non-essential requirements significantly simplify 
the discussion. More specifically, we consider a Hamilto- 
nian made of arbitrary single-body and two-body terms, 
with the interactions restricted to nearest neighbors, 
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The simulation scheme is based on an efficient decom- 
position for slightly entangled states of n-systems and on 
a protocol to efficiently update this decomposition when 
a unitary transformation is applied to one or two (nearest 
neighbor) systems, as described next. 

Efficient decomposition. — State € H-d® n can be 
expressed in terms of 0(ndX 2 ) parameters by first ex- 
panding it in a product basis, 



I*) = E " ' X Cil '- i " ® • • • ® Kn)> 



(5) 



and by then writing the d n complex coefficients Ci x ...i n in 
terms of n tensors and n — 1 vectors A^, 



ErWinWrPl'2 \[2] r [3]. 3 ...r[* / fi \ 



where each a runs from 1 to x- We refer to H for a de- 
tailed explanation on how to obtain this decomposition, 
denoted T> from now on. 

When restricted to translationally invariant states in 
one spatial dimension, decomposition T> has two main 
precedents. Except for a number of technical details, 
it corresponds to a construction introduced by Fannes, 
Nachtergaele and Werner to study the so-called finitely 
correlated states [Toj — in turn a generalization of valence 
bond states as analyzed by Affleck, Kennedy, Lieb and 
Tasaki 

EH— 

and to Ostlund and Rommer's description 
of matrix product states [12j used to analytically study 
the fixed points of the DMRG method. Here we will not 
require translational invariance and, most importantly, 
we will consider the time-dependent case. 

Vector in V contains the decreasingly ordered co- 
efficients X [ l ] > \ [ 2 ] > ■■■ > \ [ x > of the Schmidt de- 
composition of according to the splitting A : B, 
where A = [1, ■ ■ • , I], 

l*) = EAW|$L 1 -' ] }|i>L (i+1) -" ] ). (7) 

a=l 

In a generic case X grows exponentially with n. However, 
in one-dimensional settings it is sometimes possible to 
obtain a good approximation to by considering only 
the first Xe terms in Q, with X e ■C X, 
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A^H 1 -' 1 }!^-" 1 ). (8) 



This is due to the following empirical facts, involving the 
ground state \^f gr ) and low energy excitations |*fe) of 



a sufficiently regular, 2-local Hamiltonian H n Q in one 
spatial dimension. 

Observation 1: Numerical analysis shows that the 
Schmidt coefficients Aa of the ground state \^ gr ) of H n 
decay (roughly) exponentially with a, 



exp(— a) 



(9) 



Observation 2: Numerical analysis shows that during a 
low-energy evolution, as given by \*&{t)) — Ck(t)\^k)> 
the Schmidt coefficients A|* (t) also decay (roughly) expo- 
nentially with a. 

The first observation is at the root of the success of the 
DMRG in one dimension 0, 0] ■ The second observation 
implies that during the time evolution a good approxi- 
mation to |\&(f)) can be obtained by keeping only a small 
number Xe of terms in its Schmidt decomposition, leading 
to an efficient decomposition T> t . 

Simulation protocol. — Our aim is to simulate the evo- 
lution of the n systems, initially in state \^o), for a time 
T according to the Hamiltonian H n of Eq. J3J|. This 
is achieved (?) by constructing decomposition T> for state 
|*o), denoted T> , and (ii) by updating the decomposition 
T> t of the time-evolved state \^t) for increasing values of 
a discretized time t G {6, 26, ■ ■ ■ , T}, 6/T < 1. 

(i) Initialization. When |Wo) is related to the ground 
state l^gr) of Hamiltonian H n by a sufficiently local 
transformation Q |15| . 



|*o)=Q|* s 



(10) 



T>q can be obtained from decomposition V gr for [$f gr ) 
by simulating the action of Q on \^ gr ) using the present 
scheme. In turn, T> gr can be obtained through one of the 
following three methods: [if) by extracting it from the 
solution of the DMRG method (as detailed elsewhere); 
(i2) by considering any product state 

|*®} = |^M>®.~®|^M>, ($®\* gr )^0, (11) 

for which T> is always very simple, and by using the 
present scheme to simulate an evolution in imaginary 
time r according to H n , 



\Vgr) = lim 



cxp(-H n T)\$®) 
|cxp(-iJ n r)|$ 8 }| 



(12) 



or (i3) by simulating an adiabatic evolution from some 
product state to l^o) through a time-dependent 

Hamiltonian that smoothly interpolates between a local 
Hamiltonian H' n such that has as its ground state 

and Hamiltonian H n . 

(ii) Evolution. For simplicity we assume that H n does 
not depend on time [16j. Then the evolved state reads 



|* T ) =exp(-iiJ n T)|tf >. 
It is convenient to decompose H n as H n = F + G, 
F 



(13) 



X FW S X (Kf+Kf^), (14) 



even / 



G = gV] ^ ] + K 2 )> (15) 



odd I 



odd I 
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where [F®,fN] = {[G^,G^] = 0) for even (odd) I, I', 
but possibly [F, G] ^ 0. For a small S > 0, the Trotter 
expansion of order p for exp(—iH n T) reads [l7| . 



e -i(F+G)T = ^-t(F+G)S^T/S 



[f P (U FS ,U G5 )] T/s , (16) 



where 



U 



FS 



-iFS 



u. 



GS 



-iGS 



and 



fi(x,y)=xy, h{x,y) 



= x l ' 2 yx 1 ' 2 



(17) 



(18) 



for first and second order expansions (see |l7| for p = 
3,4). Eq. $nfy approximates the evolution operator 
cxp(—iH n T) by a product of 0(T/S) n-body transfor- 
mations Uf$ and Ugs, which can in turn be expressed as 
a product of two-body gates and W*p , 



u F5 = n y < 

even I 

U G s = J] W : 

odd i 



[I] 

2 ' 



[I] 

2 > 



w. 



[I] - „-iG^S 



odd I. 



(19) 
(20) 



The simulation of the time evolution (|13|l is accomplished 
by iteratively applying gates U f$ and Ugs to |^o) a num- 
ber 0(T/5) of times, and by updating decomposition T> 
at each step. If \^t) denotes the approximate evolved 
state at time t, then we have 



(21) 



where f p (UFS,UGs) consists of a product of 0(n) two- 



body gates and W^ 1 , and where we use lemma 2 in 
to update V after each of these gates. 

Errors and computational cost. - The main source of er- 
rors in the scheme are the truncation (JSJ) and the Trotter 
expansion (|16fl . We use the fidelity error 



r[l] 



e(t) = l-K* t |*t>| S 



(22) 



to measure how similar the simulated \if?t) and the exact 
|\&t) are. The truncation error e± incurred in replacing 
(J7J with |JSJ| reads 



ei 



X 

E 



[ZU2 



(23) 



On the other hand, the order-p Trotter expansion ne- 
glects corrections £2 that scale as [l8l | 



£2 



£2p T 2 



(24) 



Updating V after a two-body gate requires 0(d 3 X 3 ) basic 
operations Gates fTp* and Ugs are applied 0{T/5) 
times and each of them decomposes into about n two- 
body gates. Therefore 0(n(dx) 3 T/S) operations are re- 
quired to apply (|16[1 on \^o). If no truncation takes place, 
so that the error e is only due to the Trotter expansion, 
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FIG. 1: Propagation of a spin wave in the ferromagnetic chain 
of Eq. with n = 30 spins, B = J = 1, T = 25 and 

<5 = 0.005. The eigenvalues p a — |Aq 5 '| 2 for the reduced 
density matrix pA for half a chain are plotted as a function 
of time (thin, irregular, rapidly oscillating lines). Notice the 
exponential decay of p a in a, a seemingly common feature 
of low-energy dynamics in sufficiently regular, but otherwise 
arbitrary local models in one dimension. The time scale is 
such that the spin wave, originating at the left end of the 
open chain (see state |^o) after Eq. 12611 ') runs along the open 
chain seven times, bouncing backward each time it reaches one 
extreme spin. The error fidelity e(t) corresponds to keeping 
all 17 terms in Q and grows quadratically in the simulated 
time T (cf. Eq. 1241 1. Errors e'(t) and e"(t) correspond to 
keeping = 12 and Xe = 8 terms in Eq. (JSJ. Truncation 
errors are of the order of the neglected eigenvalues p a . 



it follows from H24[) that the number of basic operations 
or computational time T c scales as 



T c ~ n{dxf 



fl/2 



VP 



(25) 



Examples. — The numerical scheme has been tested in 
a number of one dimensional settings using matlab code 
in a pentium IV, both to (?) find an approximation to the 
ground state \^ gr ) and to (ii) simulate time evolutions of 
local perturbations of \^ gr }- The next two simple exam- 
ples aim at illustrating the performance of the method 
for order p = 2 Trotter expansions. 

(i) An approximation to the ground state of a non- 
critical spin chain, with n = 80 spins, local dimension 
d = 2, X €l = 20 and a sufficiently regular (but otherwise 
arbitrary) H n with nearest neighbor interactions can be 
obtained in half an hour by evolving some product state 
in imaginary time as in 1)12(1 , until decomposition T> T con- 
verges within an accuracy 1- |(^ r |^ r+T ,)| 2 < 10~ 10 [l9| . 

{ii) For the spin 1/2 ferromagnetic chain 



in 



1=1 



?[*+!] 



the ground state is the product state \$ gr ) 
the exact time evolution of the state |*o) 



10)* 



iiy 



(26) 



and 
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0}®" -2 can be computed efficiently, with a cost that 
grows only as 0(n 2 ), since the dynamics are confined 
to a subspace of H.2® n of dimension n(n — l)/2. In ad- 
dition, X in Eq. (JJJ is upper bounded by n/2 + 2. Thus 
we can compare the exact solution \^t) in tEH with the 
approximation \^t) obtained thorugh simulation. When 
X e = n/2 + 2, a confirmation of (|25|1 is obtained for differ- 
ent values of n, T and <5. For smaller X e , truncation errors 
of the order of the neglected Schimdt terms appear, see 
Fig. (1). 

In this paper we have shown how to efficiently sim- 
ulate low energy dynamics in one-dimensional arrays of 
n quantum systems by using 0{n) parameters. Since 
C(exp(n)) parameters are required to describe a generic 
state of n systems, this result may at first seem contradic- 
tory. However, the locality of the interactions and the ge- 
ometry of the problem make the Hamiltonian H n highly 
non-generic. It is thus conceivable that, correspondingly, 
the dynamics generated by H n are constrained to occur 
in (or can be well approximated by states from) a sub- 
manifold of TLd® n with remarkably small dimension. 

Generalizations of the present scheme to Hamiltoni- 
ans with long-range interactions, to momentum space, 
to bosonic and fermionic systems, and a specific scheme 
to address critical systems will be presented elsewhere. 
These results have also been generalized to slightly cor- 
related mixed-state dynamics, and in particular to finite 
temperature quantum many-body dynamics in one di- 
mension. 
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